module Geweke
!!DIR$ NOOPTIMIZE
	use compute_reg
	use dgp
	implicit none
	
	integer, parameter	:: npriordraws=1E5
	integer, parameter	:: nqs=19
	
	real, allocatable	:: priorstats(:,:)
	
	contains
	

	
	subroutine Gcheck
		integer		:: l,nstats,i,j,leff,ib
		real, allocatable	:: pquants(:,:,:)
		real, parameter	:: pgrid(nqs)=[(real(i)/(nqs+1),i=1,nqs)]
		real, allocatable :: chkmat(:,:), chkmat2(:,:)
		integer, parameter	:: nskip=10, nbatch=200
		real			:: cps(nbatch)
		real, allocatable	:: ix(:)
		
		call inititerate
!pop=1
Us=0
		nstats=size(tstat())
		allocate(priorstats(nstats,npriordraws),tracestats(nstats,npriordraws),pquants(nstats,nqs,2),chkmat(nstats,nqs))
		do l=1,npriordraws
			call dgp_sigma2_F
			call dgp_cindf
			call dgp_mt_F
			call dgp_F
			call dgp_s2dist
			call dgp_thdist
			call dgp_CFthdist
			call dgp_CoCthdist
			call dgp_clubid
			call dgp_CoCid
			call dgp_rhodist
			call dgp_CoCrhodist
			call dgp_cindu
			call dgp_cindCF
			call dgp_cindCoC
			call dgp_rho
			call dgp_CoCrho
			call dgp_CoCs2grid
			call dgp_clubs2grid
			call dgp_kappa2grid
			call dgp_om2
			call dgp_CoCF
			call dgp_CF
			call dgp_muU
			call dgp_U
			priorstats(:,l)=tstat()
		enddo
		call printtime
!$omp parallel do private(j)
		do i=1,nstats
			pquants(i,:,1)=quantile_v(priorstats(i,:),pgrid)
			do j=1,nqs
				pquants(i,j,2)=count(priorstats(i,:)<=pquants(i,j,1))/real(npriordraws)
			enddo
		enddo
		ix=[0.0,0.0,0.0,0.0,0.0,0.0,[(1.0+real((i-1)/2),i=1,nstats-6)]]
		print *,"quantiles of test stats"
		call mdisp(ix.clr.pquants(:,:,1))
		print *,"cdf evaluated at quantiles" 
		call mdisp(ix.clr.pquants(:,:,2))
		chkmat2=chkmat
		do l=1,npriordraws*nskip
			call draw_once
			if(mod(l,nskip)==0) then
				tracestats(:,l/nskip)=tstat()
			endif
			if(mod(l,nskip*5000)==0) then
				leff=l/nskip
!$omp parallel do private (j,cps)
				do i=1,nstats
					do j=1,nqs
						do ib=1,nbatch
							cps(ib)=count(tracestats(i,(ib-1)*leff/nbatch+1:ib*leff/nbatch)<=pquants(i,j,1))/real(ib*leff/nbatch-(ib-1)*leff/nbatch)
						enddo
						cps=cps-pquants(i,j,2)
						chkmat(i,j)=sum(cps)/sqrt(sum(cps**2))
						chkmat2(i,j)=count(tracestats(i,1:leff)<=pquants(i,j,1))/real(leff)-pquants(i,j,2)		
					enddo
				enddo
!				call mdisp(ix.clr.chkmat2)
				call mdisp(ix.clr.merge(chkmat2,0.0,abs(chkmat)>4.0))
				print *,maxval(abs(chkmat))
				print *,leff
			endif
		enddo
		print *,"Geweke check complete"
		stop
		
	end subroutine
	
	
end module